## Creation of social norms under weak institutions
#  version of March 9, 2022, last edited by florian.diekert@awi.uni-heidelberg.de
#
## Part I:
# in: csv file that has means and SE for coop rate
# out: plots of coop rate over rounds 
#
## Part II and III:
# in: csv file with all data
# out: initial beliefs plots
# out: plots of exclusion over rounds in sanctioning treatments
######-------######-------######

######-------######-------######
# clear workspace and set wd
rm(list = ls(all = TRUE))
setwd("...")
######-------######-------######

######-------######-------######
# Add libraries
library(ggplot2)
library(stargazer)
library(plm)
library(lmtest)
library(xtable)
library(readxl)
######-------######-------######



######-------######-------############-------######-------############-------######-------######
# Begin Table 3: Participants characteristics (also appendix table by treatment)
######-------######-------############-------######-------############-------######-------######

# read in data
df <- read.csv(file="all_data_R.csv",header=T)

df$treat <- NA
df$treat[df$sanc_treat==1&df$manip_treat==1] <- "hi-S"
df$treat[df$sanc_treat==0&df$manip_treat==1] <- "hi-noS"
df$treat[df$sanc_treat==0&df$manip_treat==0] <- "low-noS"
df$treat[df$sanc_treat==1&df$manip_treat==0] <- "low-S"

#x <- cbind(df$viol_area,df$viol_dagaanet,df$viol_gillnet,df$viol_license,df$viol_undersize)
#df$viol_index <- rowSums(x,na.rm=T)

varvect <- c("age","gender","earndaily","movefreq","crewsize","dagaa_gear")
treatvect <- c("hi-S","hi-noS","low-noS","low-S")

#### overview table for main text
tab <- matrix(NA,6,5)

for(i in 1:length(varvect)){
  ## define x vectors with special treatment for variables that are transformed to dummies
  x1 <- df[,which(names(df)==varvect[i])]
  if(i%in%c(3,4)){
    x1 <- ifelse(x1>1,0,1)
  }

  ## insert mean values in table
  tab[i,2] <- mean(x1,na.rm=T)
  tab[i,3] <- sd(x1,na.rm=T)
  tab[i,4] <- min(x1,na.rm=T)
  tab[i,5] <- max(x1,na.rm=T)
}
## round values
tab[1:6,2:3] <- round(tab[1:6,2:3],2)

## give names to rows and columns
tab[,1] <- c(varvect)
tab <- as.data.frame(tab)
names(tab) <- c("variable","mean","sd","min","max")

## printing of table to tex-interpretable strings
tabmat <- xtable(tab)
print(tabmat,include.rownames = F)


#### the table that splits the variables by treatment (in appendix)
tab <- matrix(NA,7,6)

for(i in 1:length(varvect)){
  ## define x vectors with special treatment for variables that are transformed to dummies
  x1 <- df[df$treat==treatvect[1],which(names(df)==varvect[i])]
  x2 <- df[df$treat==treatvect[2],which(names(df)==varvect[i])]
  x3 <- df[df$treat==treatvect[3],which(names(df)==varvect[i])]
  x4 <- df[df$treat==treatvect[4],which(names(df)==varvect[i])]
  if(i%in%c(3,4)){
    x1 <- ifelse(x1>1,0,1)
    x2 <- ifelse(x2>1,0,1)
    x3 <- ifelse(x3>1,0,1)
    x4 <- ifelse(x4>1,0,1)
  }


  ## insert mean values in table
  tab[i+1,2] <- mean(x1,na.rm=T)
  tab[i+1,3] <- mean(x2,na.rm=T)
  tab[i+1,4] <- mean(x3,na.rm=T)
  tab[i+1,5] <- mean(x4,na.rm=T)

## do the pairwise-comparisons of means, keep the min p-value
  if(i%in%c(2,3,4,6)){
    tr_1_2 <- prop.test(c(sum(x1,na.rm=T),sum(x2,na.rm=T)),c(length(x1),length(x2)))
    tr_1_3 <- prop.test(c(sum(x1,na.rm=T),sum(x3,na.rm=T)),c(length(x1),length(x3)))
    tr_1_4 <- prop.test(c(sum(x1,na.rm=T),sum(x4,na.rm=T)),c(length(x1),length(x4)))
    tr_2_3 <- prop.test(c(sum(x2,na.rm=T),sum(x3,na.rm=T)),c(length(x2),length(x3)))
    tr_2_4 <- prop.test(c(sum(x2,na.rm=T),sum(x4,na.rm=T)),c(length(x2),length(x4)))
    tr_3_4 <- prop.test(c(sum(x3,na.rm=T),sum(x4,na.rm=T)),c(length(x3),length(x4)))
  }else{
    tr_1_2 <- t.test(x1,x2)
    tr_1_3 <- t.test(x1,x3)
    tr_1_4 <- t.test(x1,x4)
    tr_2_3 <- t.test(x2,x3)
    tr_2_4 <- t.test(x2,x4)
    tr_3_4 <- t.test(x3,x4)
  }
  
  
  p.vals <- c(tr_1_2$p.value,
              tr_1_3$p.value,
              tr_1_4$p.value,
              tr_2_3$p.value,
              tr_2_4$p.value,
              tr_3_4$p.value)
  adjust.p.vals <- p.adjust(p.vals, method = "hochberg")
  minp <- min(adjust.p.vals)
  
## insert min p-value in table
tab[i+1,6] <- minp
}
## round values
tab[2:7,2:6] <- round(tab[2:7,2:6],2)



## give names to rows and columns
tab[1,] <- c("treatment","hi-S","hi-noS","low-noS","low-S","min p-val")
tab[,1] <- c("treatment",varvect)


## printing of table to tex-interpretable strings
tabmat <- xtable(tab)
print(tabmat,include.rownames = F)
######-------######-------############-------######-------############-------######-------######


######-------######-------############-------######-------############-------######-------######
# Begin Figure 2: Coop rates over rounds
######-------######-------############-------######-------############-------######-------######

# load the data
df <- read_excel("ValuesFigure_2.xlsx")

######-------######-------######
# start the figure and set graphical parameters

# start pdf or png figure. WARNING: RUN ONLY ONE OF THESE TWO LINES!
pdf("Coop_treatment_means_se.pdf",width=6,height=4.5) #width in inches

par(mar = c(4, 4, 1, 1),
    oma = c(1,0,1,1)) #bottom, left, top and right margins

plot(df$Round, df$T1m, ylim = c(.28,.63), type = "n",xaxt="n",xlab="Round",ylab="Cooperation rate")
axis(1, at=1:7, labels=c("OS","1","2","3","4","5","6"))


z <- 1 # critical value for the CI 1.96 = 95% 1.645 = 10%
x <- 2:7

## colors

alph <- 105 # alpha transperancy value (bw 0 = fully transparant and 255 = solid)

# T3: hi-exp/noS
T3.linecol <- rgb(33,113,181, maxColorValue=255, alpha = 255) # Blue
T3.bgcol <- rgb(107,174,214, maxColorValue=255, alpha = alph)
# T1: low-exp/noS
T1.linecol <- rgb(106,81,163, maxColorValue=255, alpha = 255) # Violet
T1.bgcol <- rgb(158,154,200, maxColorValue=255, alpha = alph)
# T4: hi-exp/S
T4.linecol <- rgb(35,139,69, maxColorValue=255, alpha = 255) # Green
T4.bgcol <- rgb(116,196,118, maxColorValue=255, alpha = alph)
# T2: low-exp/S
T2.linecol <- rgb(203,24,29, maxColorValue=255, alpha = 255) # Red
T2.bgcol <- rgb(251,106,74, maxColorValue=255, alpha = alph)
######-------######-------######

## make polygon where coordinates start with lower limit and then upper limit in reverse order
## make points and se for one-shot
######-------######-------######
## treatment 3 (highno sanction) in blue
L <- df$T3m[x] - z*df$T3se[x]
U <- df$T3m[x] + z*df$T3se[x]
polygon(c(x,rev(x)),c(L,rev(U)),col=T3.bgcol, border = FALSE)

segments(x0=1,y0=as.numeric(df[1,6]-df[1,7]),x1=1,y1=as.numeric(df[1,6]+df[1,7]),lwd=5,col=T3.bgcol)

######-------######-------######
## treatment 1 (low no sanction) in violet
L <- df$T1m[x] - z*df$T1se[x]
U <- df$T1m[x] + z*df$T1se[x]
polygon(c(x,rev(x)),c(L,rev(U)),col=T1.bgcol, border = FALSE)

segments(x0=1,y0=as.numeric(df[1,2]-df[1,3]),x1=1,y1=as.numeric(df[1,2]+df[1,3]),lwd=5,col=T1.bgcol)

######-------######-------######
## treatment 4 (high sanction) in green
L <- df$T4m[x] - z*df$T4se[x]
U <- df$T4m[x] + z*df$T4se[x]
polygon(c(x,rev(x)),c(L,rev(U)),col=T4.bgcol, border = FALSE)

segments(x0=1,y0=as.numeric(df[1,8]-df[1,9]),x1=1,y1=as.numeric(df[1,8]+df[1,9]),lwd=5,col=T4.bgcol)

######-------######-------######
## treatment 2 (low sanction) in red
L <- df$T2m[x] - z*df$T2se[x]
U <- df$T2m[x] + z*df$T2se[x]
polygon(c(x,rev(x)),c(L,rev(U)),col=T2.bgcol, border = FALSE)

segments(x0=1,y0=as.numeric(df[1,4]-df[1,5]),x1=1,y1=as.numeric(df[1,4]+df[1,5]),lwd=5,col=T2.bgcol)

######-------######-------######
## add lines and points
points(1,df[1,6],col=T3.linecol,pch=19, cex=1.3)
points(1,df[1,2],col=T1.linecol,pch=17, cex=1.4)
points(1,df[1,8],col=T4.linecol,pch=19, cex=1.3)
points(1,df[1,4],col=T2.linecol,pch=17, cex=1.4)

lines(df$Round[x], df$T3m[x], lwd = 3,col=T3.linecol,type = "o",pch=19,lty="dashed", cex=1.3) # blue
lines(df$Round[x], df$T1m[x], lwd = 3,col=T1.linecol,type = "o",pch=17,lty="dotted", cex=1.4) # violet
lines(df$Round[x], df$T4m[x], lwd = 3,col=T4.linecol,type = "o",pch=19, cex=1.3) # green
lines(df$Round[x], df$T2m[x], lwd = 3,col=T2.linecol,type = "o",pch=17, cex=1.4) # red

######-------######-------######
## legend all four (showing that sanctioning keeps coop at initial level: manipulation matters with sanctioning)
legend("bottomleft", legend=c("hi-S","hi-noS","low-noS","low-S"),
       col=c(T4.linecol,T3.linecol,T1.linecol,T2.linecol), lty=c(1,3,3,1),pch=c(19,19,17,17), lwd=2, cex=0.9)

######-------######-------######
## closing this figure
dev.off()

######-------######-------######
## empty workspace
rm(list = ls(all = TRUE))
######-------######-------############-------######-------######

######-------######-------############-------######-------############-------######-------######
# Figure 3 Initial beliefs
######-------######-------############-------######-------############-------######-------######

######-------######-------######
# Setting the working directory

# read in data
df <- read.csv(file="all_data_R.csv",header=T)

df$manip_treat <- ifelse(df$manip_treat==0,df$manip_treat <- "low",df$manip_treat <- "high")

######-------######-------######
## stacked barplot with all three options
# make data (need to reverse order)

df$pb_inv <- 2
df$pb_inv[df$pb==3] <- 1
df$pb_inv[df$pb==1] <- 3

df$ne_inv <- 2
df$ne_inv[df$ne==3] <- 1
df$ne_inv[df$ne==1] <- 3

table(df$ne_inv)

df$ee_inv <- 0
df$ee_inv[df$ee==0] <- 1

pb.numb <- cbind(table(df$pb_inv[df$manip_treat=="high"]),
                 table(df$pb_inv[df$manip_treat=="low"]))/length(df$pb[df$manip_treat=="low"])
ne.numb <- cbind(table(df$ne_inv[df$manip_treat=="high"]),
                 table(df$ne_inv[df$manip_treat=="low"]))/length(df$ne[df$manip_treat=="low"])
ee.numb <- cbind(table(df$ee_inv[df$manip_treat=="high"]),
                 table(df$ee_inv[df$manip_treat=="low"]))/length(df$ee[df$manip_treat=="low"])

### make plot
pdf("inibeliefs.pdf",width=7.15,height=3.2) #width in inches
par(mfrow=c(1,3),
    mar = c(4, 4, 1.5, 1)) # make the plots be closer together
# mar:  bottom, left, top, and right

barCenters <- barplot(height = pb.numb,
                      names.arg=c("high","low"),
                      #beside = true, 
                      las = 1,
                      ylim = c(0, 1.15),
                      xlab="Personal normative beliefs",
                      ylab = "Share of respondents",
                      col=c(gray(0.2),gray(0.4),gray(0.8)),
                      border = "black")


barCenters <- barplot(height = ne.numb,
                      names.arg=c("high","low"),
                      #beside = true, 
                      las = 1,
                      ylim = c(0, 1.15),
                      xlab="Normative expectations",
                      #ylab = "Share that choose 'putting points in group account' ",
                      col=c(gray(0.2),gray(0.4),gray(0.8)),
                      border = "black")


legend(mean(barCenters),1.1,fill=c(gray(0.2),gray(0.4),gray(0.8)),
       legend=c("Coop","DwoD","Defect"),horiz=TRUE,xjust=0.5,yjust=.5,
        bty="n")

barCenters <- barplot(height = ee.numb,
                      names.arg=c("high","low"),
                      #beside = true, 
                      las = 1,
                      ylim = c(0, 1.15),
                      xlab="Empirical expectations",
                      #ylab = "Share that choose 'putting points in group account' ",
                      col=c(gray(0.2),gray(0.8)),
                      border = "black")


######-------######-------######
## closing this figure
dev.off()

######-------######-------######
## empty workspace
rm(list = ls(all = TRUE))
######-------######-------############-------######-------############-------######-------######

######-------######-------############-------######-------######
# Figure 4  Evolution of empirical expectations
######-------######-------############-------######-------######

# load the data
df <- read_xlsx("ValuesFigure_2.xlsx")
df4 <- read_xlsx("ValuesFigure_4.xlsx")

df$eeT1 <- df4$T1m
df$eeT2 <- df4$T2m
df$eeT3 <- df4$T3m
df$eeT4 <- df4$T4m


######-------######-------######
# start the figure and set graphical parameters

## colors
alph <- 105 # alpha transperancy value (bw 0 = fully transparant and 255 = solid)

# T3: hi-noS
T3.linecol <- rgb(33,113,181, maxColorValue=255, alpha = 255) # Blue
# T1: low-exp/noS
T1.linecol <- rgb(106,81,163, maxColorValue=255, alpha = 255) # Violet
# T4: hi-exp/S
T4.linecol <- rgb(35,139,69, maxColorValue=255, alpha = 255) # Green
# T2: low-exp/S
T2.linecol <- rgb(203,24,29, maxColorValue=255, alpha = 255) # Red
######-------######-------######

pdf("Coop_ee_means.pdf",width=6,height=5.2) #width in inches

#pdf("Coop_ee_means_20190330.pdf",width=6,height=3.2) #width in inches, onlu hi-S and low-S

par(mfrow = c(2, 2)) # 2-by-2 grid of plots
par(oma = c(3, 4, 0, 0)) # make room (i.e. the 4's) for the overall x and y axis titles
par(mar = c(2, 2, 2, 1)) # make the plots be closer together
## bottom, left, top, and right

x <- 2:7 # which rounds to consider
ylimits <- c(.29,.62) # limits of y-axis (same in all treatments)

######-------######-------######
## erstes panel, nur y-axis ist beschriftet
## treatment 2 (private sanction) in red
plot(df$Round, df$T1m, ylim = ylimits, type = "n",xaxt="n",xlab="",ylab="",main="hi-S")
axis(1, at=1:7, labels=c("OS","1","2","3","4","5","6"))

segments(x0=-10,y0=0.3,x1=10,y1=0.3,col=gray(0.5),lwd=0.5,lty="dashed")
segments(x0=-10,y0=0.4,x1=10,y1=0.4,col=gray(0.5),lwd=0.5,lty="dashed")
segments(x0=-10,y0=0.5,x1=10,y1=0.5,col=gray(0.5),lwd=0.5,lty="dashed")
segments(x0=-10,y0=0.6,x1=10,y1=0.6,col=gray(0.5),lwd=0.5,lty="dashed")

lines(df$Round[x], df$T4m[x], lwd = 3,col=T4.linecol,type = "o",pch=19) # green
lines(df$Round[x], df$eeT4[x], lwd = 3,col=T4.linecol,type = "o",pch=19,lty="dotted") 

points(1,df$T4m[1],col=T4.linecol,pch=19)
points(1,df$eeT4[1],col=T4.linecol,pch=17)

######-------######-------######
## zweites panel, ohne jede axeis
## treatment 4 (social sanction) in green
plot(df$Round, df$T1m, ylim = ylimits, type = "n",xaxt="n",yaxt="n",xlab="",ylab="",main="low-S")
axis(1, at=1:7, labels=c("OS","1","2","3","4","5","6"))

segments(x0=-10,y0=0.3,x1=10,y1=0.3,col=gray(0.5),lwd=0.5,lty="dashed")
segments(x0=-10,y0=0.4,x1=10,y1=0.4,col=gray(0.5),lwd=0.5,lty="dashed")
segments(x0=-10,y0=0.5,x1=10,y1=0.5,col=gray(0.5),lwd=0.5,lty="dashed")
segments(x0=-10,y0=0.6,x1=10,y1=0.6,col=gray(0.5),lwd=0.5,lty="dashed")

lines(df$Round[x], df$T2m[x], lwd = 3,col=T2.linecol,type = "o",pch=17) # red
lines(df$Round[x], df$eeT2[x], lwd = 3,col=T2.linecol,type = "o",pch=17,lty="dotted") # red

points(1,df$T2m[1],col=T2.linecol,pch=19)
points(1,df$eeT2[1],col=T2.linecol,pch=17)


######-------######-------######
## drittes panel, mit allen axen (DON't run in 2x1 picture)
## treatment 1 (private no sanction) in violet
plot(df$Round, df$T1m, ylim = ylimits, type = "n",xaxt="n",xlab="",ylab="",main="hi-noS")
axis(1, at=1:7, labels=c("OS","1","2","3","4","5","6"))

segments(x0=-10,y0=0.3,x1=10,y1=0.3,col=gray(0.5),lwd=0.5,lty="dashed")
segments(x0=-10,y0=0.4,x1=10,y1=0.4,col=gray(0.5),lwd=0.5,lty="dashed")
segments(x0=-10,y0=0.5,x1=10,y1=0.5,col=gray(0.5),lwd=0.5,lty="dashed")
segments(x0=-10,y0=0.6,x1=10,y1=0.6,col=gray(0.5),lwd=0.5,lty="dashed")

lines(df$Round[x], df$eeT3[x], lwd = 3,col=T3.linecol,type = "o",pch=19,lty="dotted")  # blue
lines(df$Round[x], df$T3m[x], lwd = 3,col=T3.linecol,type = "o",pch=19,lty="solid")  # blue

points(1,df$T3m[1],col=T3.linecol,pch=19)
points(1,df$eeT3[1],col=T3.linecol,pch=17)

######-------######-------######
## viertes panel, mit nur der x-axis (DON't run in 2x1 picture)
## treatment 3 (social no sanction) in blue
plot(df$Round, df$T1m, ylim = ylimits, type = "n",xaxt="n",yaxt="n",xlab="",ylab="",main="low-noS")
axis(1, at=1:7, labels=c("OS","1","2","3","4","5","6"))

segments(x0=-10,y0=0.3,x1=10,y1=0.3,col=gray(0.5),lwd=0.5,lty="dashed")
segments(x0=-10,y0=0.4,x1=10,y1=0.4,col=gray(0.5),lwd=0.5,lty="dashed")
segments(x0=-10,y0=0.5,x1=10,y1=0.5,col=gray(0.5),lwd=0.5,lty="dashed")
segments(x0=-10,y0=0.6,x1=10,y1=0.6,col=gray(0.5),lwd=0.5,lty="dashed")

lines(df$Round[x], df$T1m[x], lwd = 3,col=T1.linecol,type = "o",pch=17,lty="solid") # violet
lines(df$Round[x], df$eeT1[x], lwd = 3,col=T1.linecol,type = "o",pch=17,lty="dotted") # violet

points(1,df$T1m[1],col=T1.linecol,pch=19)
points(1,df$eeT1[1],col=T1.linecol,pch=17)


# print the overall labels
mtext('Round', side = 1, outer = TRUE, line = 1)
mtext('Cooperation / Empirical Expect.', side = 2, outer = TRUE, line = 2)

######-------######-------######
## closing this figure
dev.off()
######-------######-------######
## empty workspace
rm(list = ls(all = TRUE))
######-------######-------############-------######-------######


######-------######-------############-------######-------######
# Figure 5 margin plots for social proximity
######-------######-------############-------######-------######

######-------######-------######

# load the data
df <- read_excel("ValuesFigure_5.xlsx")

df$affil <- c(0,0.02,0.04,0.06,0.5,0.52,0.54,0.56,1,1.02,1.04,1.06)

names(df)[names(df) == 'treatment'] <- 'treat'

z <- 1 # multiplication factor for SE
eps <- 0.008

ymin <- 0.30# min(df$coef.coop-df$SE.coop)
ymax <- 0.67#max(df$coef.coop+df$SE.coop)

treatvect <- c("hi-S","hi-noS","low-noS","low-S")

linecolvect <- c("#74C47669","#6BAED669","#9E9AC869","#FB6A4A69")
pointcolvect <- c("#238B45FF","#2171B5FF","#6A51A3FF","#CB181DFF")
linetypevect <- c(1,3,3,1)
pointtypevect <- c(19,19,17,17) 


pdf("affil_marg.pdf",width=9,height=4.5) #width in inches

par(mfrow = c(1, 2)) # 2-by-2 grid of plots
par(mar = c(4, 4, 2, 1),
    oma = c(1,0,1,1)) #bottom, left, top and right margins

## start plot cooperation

## first empty plot
plot(df$affil[df$treat=="hi-S"],df$coef.coop[df$treat=="hi-S"],
     type="n",ylim=c(ymin,ymax),xlim=c(0,1.05),
     xaxt="n",
     xlab="social proximity",ylab="predicted response",
     main="cooperation")




for(j in 1:12){
#for(j in seq(1,12,by=2)){
#for(j in c(1,3:5,7:9,11,12)){
  ## vertical lines indicating SE
  segments(x0=df$affil[j],y0=df$coef.coop[j]-z*df$SE.coop[j],
           x1=df$affil[j],y1=df$coef.coop[j]+z*df$SE.coop[j],
           col=linecolvect[which(treatvect==df$treat[j])],lwd=3)
  
  ## whiskerrs at the bottom
  segments(x0=df$affil[j]-eps,y0=df$coef.coop[j]-z*df$SE.coop[j],
           x1=df$affil[j]+eps,y1=df$coef.coop[j]-z*df$SE.coop[j],
           col=linecolvect[which(treatvect==df$treat[j])],lwd=3)
  
  ## whiskerrs at thetop
  segments(x0=df$affil[j]-eps,y0=df$coef.coop[j]+z*df$SE.coop[j],
           x1=df$affil[j]+eps,y1=df$coef.coop[j]+z*df$SE.coop[j],
           col=linecolvect[which(treatvect==df$treat[j])],lwd=3)
}


for(i in 1:4){
#for(i in 1:3){
#for(i in c(2,3)){
  lines(df$affil[df$treat==treatvect[i]],df$coef.coop[df$treat==treatvect[i]],
        type="b",lty=linetypevect[i],col=pointcolvect[i],pch=pointtypevect[i],lwd=3)
}

#legend("bottom", legend=c("hi-F","hi-noF","low-noF"),
 #      col=pointcolvect[1:3], lty=linetypevect[1:3],pch=19, lwd=2, cex=0.8)

legend("bottom", legend=c("hi-S","hi-noS","low-noS","low-S"),
       col=pointcolvect, lty=linetypevect,pch=pointtypevect, lwd=2, cex=0.75)


axis(1, at=c(0,0.5,1), labels=c("0","0.5","1"))
######-------######-------######

## start plot empirical expectations

## first empty plot
plot(df$affil[df$treat=="hi-F"],df$coef.ee[df$treat=="hi-F"],
     type="n",ylim=c(ymin,ymax),xlim=c(0,1.05),xaxt="n",
     xlab="social proximity",ylab="",
     main="empirical expectation")


j <- 1

for(j in 1:12){
#for(j in seq(1,12,by=2)){  
#for(j in c(1,3:5,7:9,11,12)){  
  ## vertical lines indicating SE
  segments(x0=df$affil[j],y0=df$coef.ee[j]-z*df$SE.ee[j],
           x1=df$affil[j],y1=df$coef.ee[j]+z*df$SE.ee[j],
           col=linecolvect[which(treatvect==df$treat[j])],lwd=3)
  
  ## whiskerrs at the bottom
  segments(x0=df$affil[j]-eps,y0=df$coef.ee[j]-z*df$SE.ee[j],
           x1=df$affil[j]+eps,y1=df$coef.ee[j]-z*df$SE.ee[j],
           col=linecolvect[which(treatvect==df$treat[j])],lwd=3)
  
  ## whiskerrs at thetop
  segments(x0=df$affil[j]-eps,y0=df$coef.ee[j]+z*df$SE.ee[j],
           x1=df$affil[j]+eps,y1=df$coef.ee[j]+z*df$SE.ee[j],
           col=linecolvect[which(treatvect==df$treat[j])],lwd=3)
}


for(i in 1:4){
#for(i in 2:3){
    lines(df$affil[df$treat==treatvect[i]],df$coef.ee[df$treat==treatvect[i]],
        type="b",lty=linetypevect[i],col=pointcolvect[i],pch=pointtypevect[i],lwd=3)
}

#legend("bottom", legend=c("hi-noF","low-noF"),
 #      col=pointcolvect[2:3], lty=linetypevect[2:3],pch=19, lwd=2, cex=0.8)

#legend("bottom", legend=c("hi-F","hi-noF","low-noF"),
  #     col=pointcolvect[1:3], lty=linetypevect[1:3],pch=19, lwd=2, cex=0.8)

legend("bottom", legend=c("hi-S","hi-noS","low-noS","low-S"),
       col=pointcolvect, lty=linetypevect,pch=pointtypevect, lwd=2, cex=0.75)

axis(1, at=c(0,0.5,1), labels=c("0","0.5","1"))


######-------######-------######
## closing this figure
dev.off()
######-------######-------######
## empty workspace
rm(list = ls(all = TRUE))
######-------######-------############-------######-------######


######-------######-------############-------######-------############-------######-------######
## The end
######-------######-------############-------######-------############-------######-------######